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INTRODUCTION:  Blast-induced  traumatic  brain  injury  (BI-TBI)  is  defined  as  the  damaging  effect  to  the  brain 
tissue  occurring  as  a  result  of  a  propagating  shock  wave  generated  by  an  explosion.  Current  TBI  models 
define  damage  in  terms  of  local  brain  accelerations  versus  time.  However,  blast  waves  can  cause  brain 
damage  by  other  mechanisms  including  excess  pressure  (leading  to  contusions),  excess  strain  (leading  to 
subdural  hematomas  and/or  diffuse  axonal  injuries),  and,  in  particular,  cavitation  effects  (leading  to  subcellular 
damage).  This  project  aims  at  the  development  of  a  realistic  and  accurate  numerical  model  that  can  capture 
the  essential  physics  of  shock  wave  propagation  and  its  interaction  with  the  skull  and  brain  tissue.  This 
information  would  enable  us  to  address  and  study  the  diverse  mechanisms  damaging  the  brain  to  better 
understand  causes  and  effects  during  traumatic  injuries. 

BODY:  Efforts  in  Project  Year  1  have  been  focused  on  several  fronts.  Our  progress  to  date  is  presented  below 
and  associated  with  each  task  (reported  here  for  reference)  in  the  approved  Statement  of  Work  . 

Task  1:  this  task  includes  the  theoretical  and  numerical  development  of  an  acoustic  and  elastic  model  to 
simulate  the  interaction  of  a  shock  wave  with  the  human  skull  and  brain.  Both  acoustic  and  elastic  wave 
propagation  components  of  the  model  are  based  on  a  system  of  first  order  partial  differential  equations 
(PDEs).  The  acoustic  model  includes  continuity  of  mass,  continuity  of  momentum,  and  a  Taylor's  series 
representation  of  the  isentropic  equation  of  state  which  links  density  and  pressure  fluctuations  and  provides 
the  source  for  the  quadratic  nonlinearity.  Tissue  specific  frequency  power-law  absorption  is  included  via  N 
independent  relaxation  mechanisms  where  the  necessary  relaxation  times  and  moduli  are  found  by  applying  a 
nonlinear  least  square  procedure.  The  elastic  model  will  be  developed  in  terms  of  coupled  isotropic  stress- 
velocity  PDEs  where  viscoelastic  processes  will  be  incorporated  via  N  relaxation  mechanisms.  Constitutive 
stress-strain  relationships  will  be  based  on  published  work.  During  this  task,  the  computational  kernel  for  both 
the  pseudospectral  method  and  the  wavelet  time  domain  method  for  adaptive  grid  refinement  based  on 
collocation  points  will  be  developed.  Specifically,  fast  and  robust  procedures  will  be  developed  for  the 
generation  of  realistic  3D  computational  domains  for  skull  and  brain;  for  the  computation  of  spatial  derivatives 
and  accurate  multi  step  time  integrators;  for  the  generation  of  multi-level  grid  adaptivity  to  dynamically  track 
the  significant  features  of  the  solution;  and  for  the  implementation  of  perfectly  matching  layers  boundary 
conditions  to  avoid  spurious  reflections  from  the  computational  boundaries.  All  subroutines  will  be  developed 
for  both  sequential  and  parallel  processing  on  scalable  single  or  multiprocessors  computer  clusters.  Parallel 
implementation  will  be  based  on  domain  decomposition  methods  using  the  Message  Passing  Interface 
protocol. 

Progress  to  date:  We  have  fully  developed  the  theoretical  model  for  both  acoustic  and  elastic  interactions  in 
heterogeneous  media  as  coupled  systems  of  first  order  partial  differential  equations.  The  model,  shown  in 
Figure  1,  is  capable  of  full  wave,  acoustic  and  elastic  propagation  in  lossy  biological  tissue.  It  accounts  for  the 
major  physical  mechanisms  that  take  place  during  the  interaction  of  shock  waves  with  tissue  and  includes 
acoustic  nonlinear  effects  to  sustain  and  propagate  the  shock  wave,  the  generation  of  both  compressional  (P) 
and  shear  (S)  waves,  and  effects  from  diffraction,  multiple  reflections,  scattering,  and  frequency  dependent 
power-law  of  absorption  which  is  characteristic  of  biological  media. 

Based  on  the  theoretical  model  we  have  also  fully  developed  the  numerical  kernels  for  nonlinear  acoustics 
and  elastic  wave  propagation  in  two-dimensional  (2D  and  2.5D  (cylindrical  symmetry))  and  fully  three- 
dimensional  (3D)  domains.  In  this  respect,  we  have  developed,  implemented,  and  debugged  source  code  in 
Fortran  90/95  for  all  the  necessary  subroutines  to  perform  the  FFT-based  pseudospectral  time  domain  method. 
The  computational  numerical  kernels  have  been  developed  in  both  sequential  and  parallel  fashion  via  the 
message  passing  interfaces  (MPI)  methods  and  domain  decomposition  techniques.  We  used  the  OpenMPI 
version  of  MPI,  which  is  a  portable,  fully  supported,  and  freely  available  implementation  of  the  MPI  standard  for 
message  passing  libraries  (http://www.open-mpi.org/).  In  particular  we  developed  algorithms  for:  1)  accurate 
pseudospectral  first-  and  second-order  spatial  derivatives  in  2D,  2.5D,  and  3D  domains  based  on  spatially 
staggered  grid.  Spatially  staggered  grid  have  been  shown  to  greatly  reduce  spurious  Fourier  oscillations  and 
Gibbs  phenomena  at  transition  boundaries  in  heterogeneous  materials,  thus  increasing  the  accuracy  and 
stability  of  the  method;  2)  a  time  staggered  fourth-order  Adam-Bashforth  method  for  time  integration.  This  is  an 
explicit  multistep  time  integration  algorithm  with  increased  stability  domain  and  accuracy,  thus  drastically 
reducing  dispersion  effects  for  a  given  time  step  compared  to  standard  methods.  It  also  allows  the  use  of 
larger  time  steps  in  the  integration  procedure  and  therefore  provides  faster  convergence  to  the  solution;  3)  the 
implementation  of  the  perfectly  matched  layers  (PML)  boundary  conditions.  PML  allow  the  acoustic/elastic 
wave  to  exit  the  computational  domain  without  numerical  spurious  reflections  at  the  boundaries  and  therefore 
greatly  aid  in  reducing  the  overall  size  of  the  computational  domain  and  the  total  number  of  computations 
without  corrupting  the  solution  with  non-physical  signals. 
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We  have  additionally  developed  a  numerical  kernel  based  on  a  wavelet-based  time  domain  method 
(WTDM)  for  dynamic  grid  adaptation  in  2D  and  2.5D  domains.  To  this  extent,  we  have  implemented 
subroutines  for:  1)  the  2D  (forward  and  inverse)  Fast  Lifted  Interpolating  Wavelet  Transform  (LIWT).  The 
single-level  fast  2D  LIWT  is  based  on  the  use  of  the  lifting  scheme  applied  to  second-generation  wavelets  and 
is  implemented  in  order  N  operations  (where  N  is  the  number  of  collocation  points  transformed)  through  the 
use  of  sequential  interpolating  filter  banks  in  a  tensor  product  approach.  M-levels  transforms  are  obtained  by 
the  recursive  application  of  the  single-level  transform;  2)  interpolation  and  spatial  partial  derivatives  routines 
based  on  cubic  splines  on  staggered  grids  where  pressure  and  velocity  components  exist  on  different, 
interlaced  spatial  and  temporal  grids;  3)  a  multi-level  grid  generation  routine  based  on  the  wavelets  coefficients 
of  the  LIWT  to  dynamically  track  the  significant  features  of  the  solution.  The  grid  generation  uses  search 
algorithms  together  with  a  thin-plate  splines  approach  to  interpolate  scattered  data  and  generate  the  adaptive, 
updated  computational  grid  at  successive  time  step  as  the  waves  propagates.  An  efficient  dynamic  grid 
solution  will  significantly  increase  the  accuracy  of  the  solution  while  drastically  reducing  both  the  computational 
time  and  resources  necessary  for  the  model. 

In  this  regard,  we  also  started  investigating  the  use  of  fast  radial  basis  functions  (RBFs)  as  an  alternative 
kernel  to  fast  Fourier  (FFTs)  and  wavelet  based  computations  for  adaptive  grid  modeling.  The  RBFs  provide 
the  advantage  of  performing  computations  on  unstructured  and  scattered  grids  and  therefore  allow  better 
handling  of  curved  interfaces  in  the  adaptive  grid  framework  in  comparison  with  the  dyadic  grid  refinement  and 
the  structured,  regular  grids  required  by  the  wavelet  and  FFTs  approaches  respectively.  The  drawback  of 
using  RBFs,  however,  lays  in  the  lack  of  very  fast  algorithms  for  non-diagonal  dense  matrix  inversion  which  are 
comparable  with  the  NIogN  operations  for  the  FFT  and  the  ill-conditioned  and  close-to-singular  nature  of  these 
dense  matrices.  We  are  currently  experimenting  with  several  approaches  to  resolve  these  drawbacks,  from 
recursive  block-decomposition  for  fast  matrix  inversion  to  several  preconditioning  methods  to  reduce 
singularity.  We  are  also  experimenting  with  the  fast  multipole  method,  a  relatively  new  approach  to  dense, 
singular  matrix  inversion  which  has  shown  promise  in  being  able  to  achieve  good  accuracy  in  close  to  order 
NIogN  computations.  We  will  thoroughly  study  this  new  approach  in  the  coming  months  to  assess  its 
performance  in  terms  of  computational  speed  and  accuracy  compared  with  the  previous  methods. 

For  the  generation  of  a  realistic  modeling  data  for  TBI  modeling  we  used  anatomical  image  slices  from  the 
Visual  Human  Project.  We  segmented,  rendered,  and  reconstructed  in  3D  the  head  section  of  the  human  man. 
Each  slice  is  segmented  according  to  the  different  anatomical  structures  present  and  each  pixel  belonging  to  a 
given  structure  is  integer-coded  (color)  in  the  database.  All  the  slices  are  then  stacked  together  to  provide  a  full 
3D  reconstruction  where  different  integers  labels  correspond  to  different  tissues/organs.  Figure  2  shows  our 
results  to  date.  We  have  the  ability  to  select  subsections  of  this  database  and  perform  computations  only  in  the 
region  of  interest.  We  are  currently  developing  subroutines  to  load  this  model  into  our  main  computation 
kernel  and  translate  the  integer-coded  tissues  to  properly  assign  acoustic  and  elastic  material  parameters. 
Furthermore,  we  are  also  in  the  process  of  developing  accurate  and  robust  interpolation  routines  to  increase 
the  resolution  of  the  3D  head  dataset  to  a  level  necessary  for  the  stability  and  accuracy  of  the  numerical 
calculations  during  shock  wave  propagation. 

Task  2:  this  task  incorporates  the  effects  of  wave  propagation  in  a  bubbly  fluid  into  the  model.  It  includes  the 
necessary  refinements  of  the  system  of  equations  that  describe  an  effective  medium  within  the  numerical 
framework  implement  in  Task  1.  As  a  shock  wave  propagates  through  a  liquid,  it  is  followed  by  a  rarefaction. 
If  the  magnitude  of  this  negative  pressure  fluctuation  exceeds  the  cavitation  threshold,  then  bubbles  will 
spontaneously  nucleate.  The  effects  of  the  bubbles  on  subsequent  shock  wave  propagation  will  be  included, 
and  the  physiological  effects  associated  with  violent  collapse  of  the  bubbles  will  be  considered.  Effective 
medium  theory  for  wave  propagation  in  a  bubbly  liquid  is  based  on  the  reduction  of  the  two  phase  media  (ie., 
gas  bubbles  and  fluid)  to  an  effective  medium  with  frequency  dependent  dispersion  and  attenuation.  The 
collection  of  bubbles  in  the  presence  of  a  sound  field  essentially  act  to  modify  the  bulk  compressibility  and 
density  of  the  fluid  to  produce  a  frequency-dependent,  complex,  bulk  compressibility  and  a  reduced  density. 
Effective  medium  parameters  change  dynamically  during  the  computations  based  on  the  bubble  sizes  and 
distributions  per  unit  volume  (the  void  fraction). 

Progress  to  date:  Towards  the  inclusion  of  cavitating  bubbles  due  to  the  shock  front,  we  have  successfully 
developed  an  analytical  representation  of  the  effective  medium  theory  for  the  different  media  present  in  the 
human  head  and  have  determined  the  necessary  effective  wave  numbers.  Effective  medium  (EM)  models 
replace  a  complex  environment  for  wave  propagation  by  an  effective  homogeneous  medium.  For  a  distribution 
of  bubbles  in  a  fluid,  an  EM  model  reduces  the  multiple  scattering  of  an  acoustic  field  to  an  estimate  of  the 
complex  wavenumber,  ke.  For  a  fluid  with  density,  sound  speed,  viscosity,  and  surface  tension  denoted  by  p,  c, 
p,  and  a,  an  EM  model  is  given  by  the  following  set  of  equations. 


5 


,2  ,2  a  2 7  a  n{ a)  da 
k~  =k+4kco2  .  v  ' - , 

o  ™\-g>  +i  ( bv+bt+ba ) 
g)q  =  |/>09? (0)—2 a/ a]  / pa2, 
bv=4  pco  I  pa2 ,  b  t=po2S(0)  I  pa~ ,  ba=co2  ka 

where  a  is  a  characteristic  radius  of  the  bubble  size  distribution  n(a),  loo  is  the  resonance  frequency  for  a 
bubble  with  radius  a,  and  bv,  bt,  and  ba  are  viscous,  thermal,  and  scattering  losses.  The  ambient  pressure  is 
po,  0  is  determined  from  the  thermal  properties  of  the  gas,  k  =  u/c  is  the  wavenumber  in  the  fluid  without 
bubbles.  To  account  for  all  orders  of  multiple  scattering,  the  scattering  loss  can  be  replaced  by  ba  =  co2kea.  For 
more  details  on  the  derivation  of  the  effective  medium  theory  please  refer  to  [1-3]. 

In  this  task,  we  are  encountering  significant  issues  in  the  integration  of  the  effective  medium  wave  numbers 
for  the  head  materials  into  the  framework  of  the  relaxation  mechanisms  approach  for  the  governing  equations. 
At  present,  our  least  square  fit  approach  for  just  two  relaxation  mechanisms,  is  not  sufficient  to  maintain 
causality  and  accuracy  in  the  attenuation  and  sound  speed  for  the  nonlinear  wave  propagation  which  results  in 
the  solution  either  showing  unphysical  behavior  or  even  growing  unstable  over  time.  We  are  currently 
investigating  these  issues  thoroughly  including  the  use  of  more  than  two  relaxation  mechanisms,  which, 
unfortunately,  will  result  in  an  increase  in  the  computational  load  and  complexity  of  the  model. 

Task  3:  this  task  includes  a  preliminary  validation  of  the  numerical  model  against  available  and  published 
benchmark  solutions  of  reduced  shock  wave  propagation  models  (such  as  the  Burgers  Equation  or  the  KZK 
model)  and  elastic  models.  Validation  can  also  be  achieved  by  comparison  with  available  experimental  data 
from  extracorporeal  shock-wave  lithotripsy.  Task  3  is  expected  to  span  a  time  period  of  about  nine  months  in 
that  any  incremental  change  to  the  mathematical  and  numerical  model  needs  to  be  revalidated. 

Progress  to  date:  We  have  validated  the  pseudospectral  and  WTDM  numerical  implementation  of  the 
individual  sequential  and  parallel  subroutines  for  spatial  derivatives  and  integration  against  an  array  of 
analytical  solutions  for  different  functional  forms  and  have  obtained  an  excellent  level  of  accuracy  which,  as 
expected,  scales  very  well  with  the  number  of  points  used  in  the  calculations.  We  have  also  validated  the  full 
sequential  and  parallel  implementation  of  the  acoustic  and  elastic  wave  model  based  on  pseudospectral 
methods  against  benchmark  analytical  solutions  for  simple  geometries  and  against  results  from  other 
published  work.  Again,  our  model  is  in  very  good  agreement  with  these  results.  Additionally,  we  also  validated 
the  nonlinear  acoustic  propagation  model  against  experimental  data  available  in  our  laboratory  through  other 
projects.  A  sample  result  showing  excellent  agreement  between  experiment  and  simulations  for  different  levels 
of  shocks  is  presented  in  Figure  3. 

Task  4:  this  task  involves  the  development  of  a  graphical  user  interface  (GUI)  to  the  numerical  kernel  to 
simplify  the  user  interaction  and  visualization  of  the  computed  results.  The  software  GUI  will  use  the  Fast 
Light  Toolkit  (FLTK).  FLTK  is  a  widgets-based,  cross-platform  C++  GUI  toolkit  which  provides  modern  GUI 
functionality  and  it  is  designed  to  be  small  and  modular  enough  to  be  statically  linked  to  other  applications. 

Progress  to  Date:  During  this  reporting  period  we  have  not  commenced  work  on  this  task. 

KEY  RESEARCH  ACCOMPLISHMENTS: 

•  We  have  fully  developed  sequential  and  parallel  kernel  versions  of  a  fully  3D  model  for  both  acoustic 
and  elastic  wave  propagation  on  fixed  computational  grids  based  on  the  pseudospectral  time  domain 
method. 

•  We  have  developed  2D  sequential  kernels  for  acoustic  wave  propagation  based  on  the  wavelet  time 
domain  method  which  include  adaptive  grid  refinement. 

•  We  have  successfully  developed  the  analytical  representation  of  the  effective  medium  theory  for  the 
different  media  present  in  the  human  head  and  determined  the  necessary  effective  wave  numbers. 

•  We  have  validated  both  the  sequential  and  parallel  implementations  of  the  acoustic  and  elastic  wave 
propagation  kernels  based  on  pseudospectral  methods  against  other  published  work. 

•  We  have  validated  the  nonlinear  acoustic  propagation  model  against  experimental  data  available  in  our 
laboratory  through  other  projects. 

•  We  have  completed  segmentation,  labeling,  and  reconstruction  of  the  head  region  of  the  Visual  Human 
Project  man. 
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REPORTABLE  OUTCOMES: 


Poster  presentation:  “Comprehensive  Three  Dimensional  Model  of  Shock  Wave-Brain  Interactions  in  Blast 
Induced  Traumatic  Brain  Injury”,  Department  of  Defense  Third  Military  Health  Research  Forum,  Hallmark 
Crown  Center,  Kansas  City,  Missouri,  August  31  -  September  3,  2009 

CONCLUSION:  We  have  developed  and  validated  a  fully  3D  acoustic  and  elastic  wave  propagation  model  to 
investigate  the  effects  of  blast-induced  shock  waves  to  the  head  and  brain  regions.  The  modeling  domain  is 
based  on  the  segmentation,  labeling,  and  3D  reconstruction  of  real  anatomical  slices  from  the  Visual  Human 
Project.  Inclusion  of  cavitating  bubble  effects,  extension  of  the  numerical  algorithms  to  adaptive  grid 
refinement  during  the  calculations,  and  the  addition  of  a  graphical  user  interface  to  the  software  are  currently  in 
progress  and  under  development. 
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SUPPORTING  DATA: 


Full-wave,  nonlinear  acoustic  propagation  in  lossy  tissue 
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Figure  1  -  Comprehensive  coupled  PDE  systems  for  modeling  shock  wave  propagation  in  biological  media.  A.  Full 
wave,  acoustic  propagation  in  lossy  tissue  accounting  for  nonlinearities,  diffraction,  multiple  reflections,  scattering,  and 
frequency  dependent  power-law  of  absorption  through  N  relaxation  mechanisms.  B.  Comprehensive  system  of 
equations  for  viscoelastic  propagation  (stress-velocity  formulation).  Full  wave  propagation  in  lossy  solids  accounting 
for  compressional  (P)  and  shear  (S)  waves,  diffraction,  multiple  reflections,  scattering,  and  frequency  dependent 
power-law  of  absorption  through  N  independent  relaxation  mechanisms.  Equations  (1)-(2)  and  Equations  (3)-(4)  refer 
to  the  dynamics  of  the  diagonal  (i=j)  and  off  diagonal  (i^j)  components  of  the  stress  tensor  and  relaxation  memory 
variables,  respectively.  Equation  (5)  governs  the  conservation  of  momentum  and  completes  the  system  of  first-order 
PDEs.  List  of  symbols  (A):  p  acoustic  pressure,  v  acoustic  particle  velocity,  p0  material  density,  k  compressibility 
moduli,  t  relaxation  times,  B/A  nonlinear  coefficient,  S  state  variables.  List  of  symbols  (B):  cr(/  is  the  ijth  component  of 

the  symmetric  stress  tensor;  v.  denotes  the  ith  component  of  the  velocity;  r  ~  are  the  memory  variables;  rf  and  xs£ 
are  the  strain  relaxation  times  for  the  P-waves  and  the  S-waves,  respectively;  t a  is  the  stress  relaxation  time  for  both 
P-  and  S-waves;  //  is  the  relaxation  modulus  corresponding  to  S-waves  and  it  is  the  analog  of  the  Lame  constant  p  in 
the  classical  elastic  case,  k  is  the  relaxation  modulus  corresponding  to  P-waves  analogous  to  A+2\i  in  the  classical 
elastic  case;  p  is  the  density. 
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Figure  2.  3D  segmentation  and  reconstruction  of  the  head  region  of  the  Visual  Human  Project  (VHP).  Left: 
Segmentation  and  color-coded  fill  of  the  major  tissue/structures  in  the  original  anatomical  slice  of  the  VHP  as 
provided  by  the  NIH  database.  Center:  3D  contours  stack  of  all  the  segmented  slices.  Right:  Full  3D  rendering  of 
the  head  section  showing  the  different  structures  and  the  ability  to  chose  the  intersecting  planes  of  the  interested 
volume. 
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Nonlinear  waveforms  modeled  < - )  and  measured  ( - )  at  focus  at  different  power  levels 

1  mm  behind  an  absorptive  (0.7  dB/cm/MHz)  phantom  of  31  mm  thickness 


Figure  3.  Validation  of  the  acoustic  modeling  versus  experiment  in  absorptive  phantom.  Nonlinear  waveforms  with 
shocks  were  measured  and  modeled  at  focus  in  less  than  0.5  mm  beyond  a  phantom  of  31  mm  thickness  and 
found  in  excellent  agreement. _ 
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